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Abstract 

We study scaling behavior and phase diagram of a 'P7~-symmetric non-Hermitian Bose-Hubbard model. In the free 
interaction case, using both analytical and numerical approaches, the metric operator for many-particle is constructed. 
The derived properties of the metric operator, similarity matrix and equivalent Hamiltonian reflect the fact that all the 
matrix elements change dramatically with diverging derivatives near the exceptional point. In the nonzero interaction 
case, it is found that even small on-site interaction can break the TT symmetry drastically. It is demonstrated that 
the scaling law can be established for the exceptional point in both small and large interaction limit. Based on 
perturbation and numerical methods, we also find that the phase diagram shows rich structure: there exist multiple 
regions of unbroken VT symmetry. 
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1. Introduction 

Non-Hermitian Hamiltonian is traditionally used to describe open system phenomenologically. It has profound 
applications in nuclear physics, quantum transport, quantum chemistry, as well as in quantum optics yj]. Since 
the discovery of a parity-time (PT~) symmetric non-Hermitian Hamiltonian can still have an entirely real spectrum 



extensive efforts were paid to the pseudo-Hermitian quantum theory fl2l |3U4ll5Ll6ll7Ll8ll9ll !CtlllllT^n3,[T^ITHfT6l[T7 . 
18, 2^, 21, 22, 23, 24, 25], which paved the way to our understanding of the connection between non-Hermitian 
systems and the real physical world. In general, a PT -symmetric non-Hermitian Hamiltonian has unbroken as well 
as broken VT -symmetric phases, the phase boundary is referred to as the exceptional points (EPs). Studies of the EPs 
were presented theoretically and experimentally over a decade ago J3.S0]. Recently, the experimental realization 
of !PT~-symmetric systems in optics were suggested through creating a medium with alternating regions of gain and 
loss |2^, 21, 22], in which the complex refractive index satisfies the condition V (x) — V* (-x) and PT symmetry 
breaking was observed 1123 1. 

One of the characteristic features of the !P7~-symmetric system is the ubiquitous phase diagram which depicts the 
symmetry of the eigenfunctions and the reality of the spectrum ll ill . The phase separation arises from the fact that 
although H and the TT operator commute, the eigenstates of H may or may not be eigenstates of the VT operator, 
since the PT operator is not linear. In the broken !PT-symmetric phase the spectrum becomes partially or completely 
complex, while in the unbroken !P7~-symmetric phase both H and VT share the same set of eigenvectors and the 
spectrum is entirely real. Recently, the phase diagram of a lattice model has been investigated. It is shown that the 
critical point is sensitive to the distribution of the coupling constant and on-site potential (2' 



27, 28, 2 



In this paper, we investigate the effect of on-site interaction on the phase boundary of a !PT-symmetric Bose- 
Hubbard system. Our approach is based on our previous work in Ref. 12411 . where we have systematically investigated 
an Af-site tight-binding chain with a pair of conjugate imaginary potentials +iy located at edges. Here we will gen- 
eralize this description by considering many-particle system and adding the on-site Hubbard interaction U. In the 
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free interaction case, many-particle eigenstates are obtained in aid of the single-particle solutions. We also construct 
the metric operator to investigate the Hermitian counterpart and observables in the framework of complex quantum 
mechanics. In nonzero U case, we restrict our attention to the influence of the nonlinear on-site interaction U on the 
boundary between unbroken and broken !PT-symmetric phases. Exact Bethe ansatz solution and numerical results 
show that small on-site interaction can reduce the critical point y c drastically. Moreover, numerical results show that 
there exist multiple regions of unbroken VT symmetry and the number of such regions increases as the system size 
N increases. 

This paper is organized as follows. Section|2]describes the Hamiltonian of a VT -symmetric Bose-Hubbard model. 
In Section [3] we focus on the interaction-free case. Based on the single-particle solutions, we construct the many- 
particle eigenstates and metric operator to study the Hermitian counterpart and observables. Section |4]is devoted to 
the case of nonzero interaction. Based on the approximation solutions, we investigate the critical scaling behavior 
and the phase diagram. Our findings are briefly summarized and the physical relevance of the model and results are 
discussed in Section 



2. !PT-Symmetric Bose-Hubbard model 



The Bose-Hubbard model gives an approximate description of the physics of interacting bosons on a lattice. Since 
it embodies essential features of ultracold atoms in optical lattices, the Bose-Hubbard model plays an important role 
in quantum many -body physics 1 3(1 31 ]. Optical realization of two-site Bose-Hubbard model in coupled cavity arrays 
and waveguides have been proposed D2[ l33l l34ll . It is worth notin g th at non-Hermitian Bose-Hubbard dimer has 

' \. Theoretical investigations on two site 



36, 37 



attracted enormous research attention in recent years []35 
open Bose-Hubbard system was firstly presented in 0511 . For a VT -symmetric non-Hermitian Bose-Hubbard dimer, 
the spectrum and the exceptional points were studied in 113611 . After that, dynamics in a leaking double well trap 
described by non-Hermitian Bose-Hubbard Hamiltonian with additional decay term was investigated under the mean 
field approximation [371- Through dynamical study of Bose-Einstein condensed gases, it was shown that imaginary 
periodic potential may induce perfect quantum coherence between two different condensates 1I39I1 . The realization of 
such open system can be put into practice as a BEC in a double well trap, where the condensate could escape from the 
traps via tunneling. Most investigations mainly focus on the Bose-Hubbard model with effective decay term in one 
site. However, it should be noticed that non-Hermitian Bose-Hubbard dimer with complex coupling terms has also at- 
tracted some attention, the decay of quantum states could be controlled by modulating the particle-particle interaction 
strength and the dissipation in the tunneling process f[4Qh . On the other hand, the Bose-Hubbard model with particle 
loss was investigated in an alternative way through employing Lindblad master equation [41], which phenomeno- 
logically describes non-unitary evolution of an open system 114211 . Recently, !P7~-symmetric quantum Liouvillean 
dynamics is also investigated 14311 . In this paper, we investigate the property of a non-Hermitian Hamiltonian in the 
framework of quantum mechanics. 

Nevertheless, although there have been no experiments to show clearly and definitively that a finite non-Hermitian 
Hamiltonian do exist in nature, many interesting features have been observed in non-Hermitian optical systems, such 
as double refraction, power oscillations, nonreciprocal phenomenon, etc. |21|,|22|]. So far, most contributions 

to pseudo-Hermitian quantum theory were for the single particle problem. Particularly, a two-mode P F-sy mmetric 
non-Hermitian Bose-Hubbard system with an imaginary potential on the edge has been investigated 0611 . In this 
paper, we focus on the influence of on-site Hubbard interaction, not restricted to a dimer but to two-particle problem 
of an A^-site Bose-Hubbard system. We mainly study the !P7~-symmetric non-Hermitian Bose-Hubbard system. The 
Hamiltonian reads 

H = -J 2^ ( a ]«/+i +H.c.) + jYjrfrf + iy( - ni ~~ ^ 
1=1 1=1 

where a\ is the creation operator of the boson at the Zth site and the tunneling strength is denoted by J. The on-site 
interaction strength and the on-site potential are denoted by U and iy, respectively. H is a VT -symmetric Hamiltonian, 
i.e., [VT, H] = 0, where the action of the parity operator V is defined by V : I — > N + 1 - / and the time-reversal 
operator T by T : i — > —i. Both single-particle solution and the critical point y c for interaction-free Hamiltonian 
jjhee _ jj^jj _ q-j naye been obtained explicitly in our previous study Il24ll . The main goal of the present work is to 
study the influence of the nonlinear interaction on the features of the system. 
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3. Interaction-free system 



When dealing with a PT -symmetric non-Hermitian system, to our knowledge, most researchers concerned about 
the single-particle problem, since it is believed that the extension of this study to a many-body problem is straight- 
forward. Nevertheless, due to the particular formalism of non-Hermitian quantum mechanics, it is worthwhile to 
investigate the many-particle system with U = 0. In the following, we will extend the obtained results for single 
particle to the case of many-particle system with zero U. 

3.1. Many-particle solutions 

In a non-Hermitian system, although the particle probability is no long conservative, the particle number 



N„ 



a. 



/=! 



still shares the common eigenfunctions with the Hamiltonian due to the commutation relation 



[N p ,H bee ] = [N p ,H] = 0. 



(2) 



(3) 



This fact indicates that the proper inner product should accord with the conservation of particle number. Therefore the 
eigenstates of H tlee or H can be obtained in each invariant subspace V Np , which is spanned by the occupation number 
basis 



\ni,n 2 ,...,n N ) = ]~ [|«,>, 



(4) 



with N p = n h where \n{) = (a) )"'/ Vulvae). Notice that [\n\,ri2, «jv)} is orthonormal set under the Dirac inner 
product. According to our previous work [01, the single-particle solutions {|<^+^} and {|<^.\}, which are eigenfunctions 
of the systems H flee and (H free ) T , i.e. \(p+J = Yii fl a ] l vac X \<P-) - Higl®] |vac). They can construct the biorthogonal 



basis set, i.e., 



2 ■/*<&* 



(5) 



here f[ = cf> k + (I), g' k = <p k _ (I) have the form 



e 



m-No) _ 



-W+No) 



where 



J[l + l^±| 2 ] sin (Nk) /sink- 2N£ ± 
ye' k + iJ 



ye 



+ U 



(6) 



(7) 



and iVo = (N + 1) /2. The symmetries of the wavefunctions and the spectrum reveal that there are two phases, 
unbroken and broken phase, which are separated by the critical point y c , 



Jc = 



±J, N = 2n 

+J J*2±, N = In + 1 



where n = 1,2, ... In the unbroken region with \y\ < \y c \, all the solutions possess !P7~ symmetry 

prfi={fr i - i y=fL 



(8) 



(9) 
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and the spectrum is entirely real. In the following, we will demonstrate that the above analysis can be extended to 
many-particle sector. Actually, one can define the operators in k space in the form of 



d k 



(10) 



a k = lj\g k ) a u 

i 



(11) 



which obey the standard bosonic commutation relations 

[at, at'] = S kk >, 
[a k ,a k >] = [a k ,a k >] = 0. 

Then the Hamiltonian H free can be written as the diagonal form 

HfKe =Yi £k ~ akClk ' (12) 

where e k = -27 cos k is real. With respect to the canonical commutation relations of ( fTTI ). the Hamiltonian in the 
form of (O can be regarded as the term of the so-called second quantization representation. Defining the occupation 
number state in £-space as 

K) = , |vac) , 



jn ki \ 



(«[)"'■ 
K) = —j= I vac), 



(13) 



which satisfy 



a kl m, 



) = V"A, + 1 n ki + 1), 



a k \n k ) = y/n^n^ - 1), 
a\, \n k ) = yjn ki + 1 \n k . + 1) . 
a\, \n k .) = V«^| W A', - 1) • 
Then, the eigenstates in the subspace V Np of H bee and (H hee )^ read 

N 



(14) 



\n ki ,n kl , ...,n kN ) = ]~~[ |n A . ), 
i=i 

A? 

\n kl ,n k2 , ...,n kN ) = Y\\ n k)> 



(15) 



respectively. They correspond to the same eigenvale as 

N 



E(n kl ,n k2 ,...,n kN ) = ^n kl e kl . (16) 



and the total particle number as N p = 2/=i n k r Equivalently, we have 



H \n kl ,n kl , ...,n kli ) = E (n kl ,n kl , ...,n kN ) \n kl ,n k2 , ...,n kfi ), (17) 
(// tree )' \n kx ,n k2 ,...,n kN ) = E(n kl ,n k2 , ...,n kN ) \n kl ,n k2 , ...,n kli ) . (18) 



Thus we conclude that for many -particle case, the phase boundary is still at y c . Notice that the eigenstates {\n kl , n kl , n kti }}, 
n kl ,n kl , n kN )} construct a biorthogonal set instead of the set {|«i,«2> n N)} under the Dirac inner product. 
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: y, and r = yjj 2 - y- 12. These matrices satisfy the relations \22\ . The 



Table 1 : The matrix representation of the metric operator rj, similarity matrix p, and the equivalent Hermitian Hamiltonian h for systems with 
N = 2, 3, and 4 are listed. Here we denote A = Jfi 



■ J, «± = 

analytical expression for the cases of N = 2, 3 and the numerical plot in Fig. ^for the case of N = 4 show that the derivatives of all the matrix 
elements with respect to y diverge at the exceptional points. 
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3.2. Metric and Hermitian counterpart 

According to the quasi-Hermitian quantum mechanics, a bounded positive-definite Hermitian operator n in each 
invariant subspace can be constructed |4|] via the eigenstates of (// free )^ as 

r] = ^ \n kl , n kl , ...,n kfl ) (n kl ,n k2 , ...,n kfl \ , (19) 

which is called the metric operator to define the biorthogonal inner product. The 77-metric operator inner product leads 
to a unitary time evolution Ji, H]. Here {n ki } denotes all the possible states with 2*. n k, — N p . One can see that the 



metric operator fulfils 

tjH^t]- 1 = (H fKs y, (20) 

and thus can be employed to construct a Hermitian Hamiltonian h that possesses the same spectrum as // flee . Actually, 
the matrix representation of 77 based on the orthonormal basis under the Dirac inner product, says shows that it is 
a Hermitian matrix. Furthermore, let p = ^jrj be the unique positive-definite square root of 77. Then the Hermitian op- 
erator p acts as a similarity transformation to map the non-Hermitian Hamiltonian H flee onto its equivalent Hermitian 
counterpart h by 

h=pHp~\ (21) 

To demonstrate such a procedure we take the small size systems as examples. In the following, we consider the 
Hamiltonian matrices in single-particle subspace for chain systems with N = 2, 3, and 4. We derive the explicit 
forms of metric operator 77, similarity matrix p, and Hermitian counterpart h for non-Hermitian Hamiltonian H^. The 
matrices 77, p, and h for systems = 2, 3 are expressed in analytical forms in Table Q] while the ones for = 4 are 
plotted in Fig. Q] It is noticed that they satisfy the following relations 

<RifR = rf=rf x , (22a) 
<Rp<R = p*=p~\ (22b) 
PTrfPT = T],PTpPT = p, (22c) 
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Figure 1 : Plots of the derivatives with respect to y of matrix elements of r\ (left panel), p (middle panel), and /; (right panel) for the case of N = 4 
near the exceptional point y — > y c = 1 (taking J = 1), The numerical results are obtained by exact diagonalization, where (;',/) denotes the index 
of the matrix elements. It shows that the derivatives of the elements diverge as y — > y c . 



where matrix 7? is defined as % (in, n) = (-1)"' S mn . 

All the matrices have the common features: the derivatives of them with respect to y diverge at the exceptional 
points. This result is not surprising since there is at least a pair of energy levels exhibit repulsion characteristic near 
y c . Nevertheless, we notice that the derivative of the original non-Hermitian Hamiltonian H free is always finite, which 
reveals the essential difference between a pseudo-Hermitian Hamiltonian and its equivalent Hermitian counterpart. 
The physics of h near the exceptional point is also obvious: the coupling constants of Hermitian counterpart change 
dramatically with diverging derivatives. 

3.3. Observables 

Another theoretical interest in non-Hermitian VT -symmetric system is that the unitary evolution can be obtained 
by introducing metric operator. In this section, we will illustrate the basic ideas via the above analytical solution. 

By introducing ^-metric operator inner product, = (-| tjV), time evolution can be expressed in a unitary way 
and also a fully consistent quantum theory can be established JHQj]]. Accordingly, the physical observables O with 
respect to the metric operator 77 can be constructed to meet the relation 

r,Ori- { = O 1 . (23) 

We examine the total particle number operator. It is defined as 

N 

1=1 k, 

In the invariant subspace V Np , we have [N p , H free ] = 0, which allows the eigen equations of the operators in the form 

of 

N p \n h ,n kl ,...,n kN ) = N p \n kl ,n h , ...,n kN ), 

N p \n kl , n kl , n klr ) = N p \n ki , n kl , n klf ) , 

and 

N p \n x ,n 2 ,...,n N ) = N p \n l3 n % , ...,n N ) , (26) 
These indicate that the Hermitian operator N p is an observable. Alternatively, this can be proved in the framework of 
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non-Hermitian quantum mechanics. Actually, we have 



TftpV ' = X \ n k l ,nk 2 ,--,n k J{n kl ,n kl ,...,n kN \ (27) 

X 2 ' ' " Uk 'N^ nk \ ' " 4 2' -' 

{n k ,} 

= n p X l"*i' n *2' •••>«*«) ( n ki,n k2 , -,n kN \ 



= N 1 



since the biorthogonal basis satisfies 



z 



n kl ,n kl , ...,n kfi ) (n kl ,n k2 , ...,n kll \ = 1. (28) 



K-) 

Accordingly, for operator h kl = a kl a kn we also have 

Wk,ri~ l =^\n kl ,n k2 ,...,n kN ){n k] ,n k2 ,...,n kN \ (29) 

K) 



=4' 

in order to obtain d27} and (f29) we used (fl?t . (l25l l, and (l28l . 

Then we conclude that the two types of particle number operators, N p = Yif=\ a \ a i an d n kl = a kl a kl , are both 
observables. Besides, the Hamiltonian H flee itself and the metric operator rj are also observables. Here, we would 
like to clarify that N p and 77 are both Hermitian operators, but H hee and n k] are non-Hermitian operators. However, 
operators a\ a, and aj, a k] are no longer observables 13L which are proved through a simple illustration in the 



Appendix A 



Here we want to stress that, the term "observable" is specific to the non-Hermitian quantum mechanics framework, 
which differs from that in traditional quantum mechanics. It is still controversial for the interpretation of the observ- 
able. As pointed above, N p is a good quantum number, or say, the obtained eigenstates of H are also the eigenstates 
of the total particle number N p . This guarantees a unitary time evolution if the 77-metric operator inner product is 
taken. However, the Dirac expectation value of the particle number is not conservative under time evolution, which 
seems to be expected. This is basically caused by the fact that the eigenstates of !PT-symmetric Hamiltonian are non- 
orthogonal under the Dirac inner product. On the other hand, all Hermitian operators are observables in Hermitian 
quantum physics. For example, operator a! a,- is an observable in Hermitian quantum mechanics but not regarded as 
an observable according to the non-Hermitian theory. Nevertheless, it is worth to mention that the observation of the 



non-Hermitian behavior in experiment, e.g., the power oscillation phenomenon 12 111 02311 . is based on the distribution 
of Dirac expectation value for aja*. 



4. Nonzero interaction system 



Now we turn to investigate scaling behavior and phase diagram of the system at nonzero U. The boundary of the 
phase is the main character for a non-Hermitian system. So far most studies dealt with the noninteracting system. 
For a non-Hermitian lattice model, it is shown that the critical point is sensitive to the distribution of the coupling 
constant and on-site potential 11261127112811 . It indicates that the inhomogeneity of a noninteracting system may shrink 
the unbroken region of VT symmetry. From the point of view of mean field theory, on-site interaction takes the role 



7 



of on-site potentials in some sense. Thus it is presumable that a nonzero U may shift the critical point. In most cases 
of nonzero U, an exact solution is hard to obtain. In this paper, we only consider the two-particle case within some 
specific parameter areas. 

4.1. Solutions for nonzero U 

The two-particle Bethe ansatz solution 

k*i,fe> = ^]fkiM (h,h)a\a] 2 |vac> (30) 

where the explicit form of f^fo (h,h) can be expressed as 

fhteihJz) = A(k l ,k 2 )e ik,l,+ik2 ' 2 +A(k 2 ,ki)e ik2h+iklh + A (Jti, -Jk 2 ) c* : ' 1_<fefc +A{k 2 ,-ki)e ikh ~ iklh (31) 
+A (-k u k 2 )e~ iklh+ik2h + A (-k 2 , ki)e~ ik * h+iklh + A (-k u -k 2 ) e^ Clh ^ h + A (-k 2 , -k x ) 

Based on the stationary Schrodinger equation 

H\li kiM ) = E{k u k 2 )\li kiM ) , (32) 

quasimomenta k\ and k 2 satisfy the equations 

(/2 + y V'2M e -WD Q Q 

= — (33) 

(j2 + y 2 e -i2k l J e ik l (N + l) G—G+- K ' 

k x ^k 2 . (34) 

where 

Go-, T > = 2Ji sin k\ + cr2Ji sin k 2 + cr'U (35) 

with cr, cr' = +. Hereafter k\ <-> k 2 denotes the corresponding equation by exchanging k\ and k 2 . The quasimomenta 
k\, k 2 and amplitudes A {ki,k 2 ) can be determined by d32l and the proper definition of inner product according to 
the VT -symmetric quantum theory. The corresponding eigenvalue is E{k\ , k 2 ) — -2 J (cos k\ + cos k 2 ) , the reality of 
E(k\,k 2 ) determines the phase diagram of the system. It is obvious that d33l and d34l l are invariant under k\ — > k 2 , 
k 2 — > k\\ and also under k\ — > -ij, A:2 — > -A^. We rewritten (l33l and (l34l explicitly as 

{VI J) sin Jti {cos [Jti (AT + 1)] + (r/7) 2 cos \k\ (N - l)]j 

+ {sin[fc 1 (A^+ 1)] + (y/7) 2 sin (A^ - 1)]} , (36) 
x [sin 2 k 2 - sin 2 h + (U/2J) 2 ] = 

k x «-» fe. (37) 

Although the analytical solutions of (l36l l and (|37] | are hard to obtain, approximate solutions within certain ranges of 
the parameters y, J, and U may shed light on the influence of U on the phase boundary. 

4.2. Solutions at the point y — J 

We start with the solution of an even system at the point y — J, which is the exceptional point for the system of 
U = 0. The influence of on-site interaction on the phase boundary can be qualitatively revealed. Taking y = J, (f36b 
and (l37l i are reduced to 

t/ sin cos (MO + J [(U/2J) 2 + cos 2 /fc,] x sinfeAO = 0,k 2 = tt/2, (38) 

h «-» k 2 . (39) 

respectively. We notice that k\ = n/2 is the solution of (l38l for U — 0. Then for small J/ case, the solutions should 
have the form (tt/2,jt/2 + 8 2 ) and (n/2 + 0\,n/2). For sufficient small U, taking the approximations sin^i^AO a 61,2^ 
and cos (9i2N) ~ 1, the critical equation reduces to 

6 3 l2 + (U/2J) 2 u2 + U/(JN) = (40) 
8 



which has one real root and two non-real complex conjugate roots, since the discriminant of the cubic equation 
A = [U/ (2JN)] 2 +[(U/2J) 2 /3] 3 > 0. Furthermore, one can get the solution of the cubic equation by routine. 
In order to obtain a concise expression of 0\2, we simply ignore the term of U 2 in the cubic equation, and then 
obtain 6^2 = -^[UjJN, (1 + i yf§)i/U/8JN. The corresponding complex conjugate eigenvalues are E± = 2/ sin #1,2 
~ (1 + i^)y/UJ 2 /N. Then we can conclude that point y = J is in the broken VT -symmetric region in the presence 
of on-site interaction, which shrinks the unbroken region of VT symmetry. Note that for a given N, the eigenvalues 
E ± become further away from real values as U grows. This agrees with the numerical simulation for phase diagram 
of the finite size systems. 

4.3. Exceptional points and scaling behavior 

Now we focus on the phase boundary of the system with small U. It is presumable that one pair of coalescing 
eigenstates have the quasimomenta k\ 2 = x/2 + <5i,2 with ]<5i 7 2j 1- The original critical ((36) and (l37t are reduced to 

(U/J) cos 81 {sin [5 1 (N + 1)] - (y/J) 2 sin [81 (N - 1)]} 

- {cos (N + 1)] - (y/J) 2 cos (N - 1)]} , (41) 
x [cos 2 5 2 - cos 2 61 + (U/2J) 2 ] = 

(5 2 *->c5i. (42) 

Furthermore, under the approximation j <5 1 ,2 j N <sz I and ignoring the term of U 2 , (f36l > and ( l37l l are reduced to polyno- 
mial equations 

(6 2 1 -8l)(t-6f)-(N 2 Z+l)6 l u = 0, (43) 
8 2 <r*8 u (44) 

where we have defined 

J 2 - y 2 U 
J 1 + y l J 

Eliminating 82, one can obtain the equation about 8\ in the form of 

f(Sx)=8\- (8\ + u {(N 2 + l)6{- ?8\ + e = 0, (46) 



which solution determines the eigenvalues. As pointed out in Ref. II24H . when the eigenstates turn to coalescence at the 
critical point y c , f(8\) should also satisfy the equation df(5\) /d8i = 0. Eliminating 8\ from (l46l i andd/(5i) /d<5i = 0, 
we have 

3 3 u 2 ({N 2 + l) -2 8 ^ = (47) 
under the condition \(\ N 2 <«c 1 . Then we can obtain y c approximately as 



yc ^ J ^Y^" H1 ~^' (48) 

where f3 = (3/8) -\j2NU 2 1 J 2 . It shows that the exceptional point exhibits an interaction sensitivity and undergoes 
dramatic changes following the change of the Hubbard interaction. Substituting y c into d43l and (PPfl) . we have 



1 ±/V2 ,I~U 

01,2 = 

which leads to the critical eigenvalue 



E c = 2 J (sin 81 + sin 8 2) ~ 2J ^]jJ]j ■ (50) 




Figure 2: (Color online) Plots of i5 1 1 . i52\ , i53\ and the corresponding numerical simulation obtained by exact diagonalization. In (a) and (b), the 
blue crosses, circles and squares indicate the numerical results for the cases of U = 10~ 4 , 10~ 5 and 10~ 6 , respectively. The black lines are the plots 
of the corresponding analytical expressions. It shows that they are in agreement with each other for U = 10~ 5 and 10~ 6 . In the case of U = 10~ 4 , a 
slight deviation appears for large N. From (c), it is observed that the numerical result accords with the analytical expression very well for N = 10 
and 20. A slight deviation appears for large U in the N = 40 system. 
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Figure 3: Phase diagrams obtained by exact diagonalization for finite systems with N = 10, 20, and 40. The shadow area indicates the region 
where the spectrum is entirely real. It shows that as N increases the unbroken TT -symmetric region becomes narrow. There are multiple unbroken 
PT -symmetric regions appear as N increases. 



At U — 0, d48l l and d50b reproduce the obtained results in our previous work [24]: y c = J and E c = 0, respectively 



Furthermore, it is observed that for small U and finite N, the critical quantities y c and E c can be expressed as 

1. . , 3 /CA 2/3 



3 lUv 1 

ln(l- 7c /7) * -lnjV + ln[— (y) ], (51) 
\n(E c /J) * _-lnjV + ln[2 2/3 ( 7 j ], (52) 



which shows the similar dependence on the size of the system. According to the finite-size scaling ansatz 114411 . the 



critical behavior can be extracted from the above-mentioned finite samples. Then combining ( |48| ) with ( l50l l leads to 

J 2 -r 2 3 

Ec-fr^ = 7 U, (53) 
J 2 + H 4 

which is a universal scaling law for such a phase transition in small U limit. To verify and demonstrate the above 
analysis, numerical simulations are performed to investigate the scaling behavior. We compute the quantities y c and 
E c for finite systems, which are plotted in Fig. [2] as comparison with the analytical results dBTT l. d52l , and (l53l . It 
shows that for small U, they are in agreement with each other. It is worthy to note that the analytical expressions in 
( 1511 ), ( 1521 , and (|53]l are obtained under the condition UN 2 /J «c 2 4 / V3 1 (obtained from |^ lj2 | N «c 1 and |^| N 2 <sc 1). 



10 




Figure 4: Plots of function Z (K), insert is Z (K) in the region (0, n). (a) N = 10, (b) N = 20, (c) N = 40, (d) W = 200, (e) N = 1000, (f)N = 5000. 
The number of roots of Z(K) = is 1, 3, 5, 1 1, 27, 63, the number of unbroken TT symmetric regions is 1,2, 3, 6, 14, 32, respectively. The black 
line, Z(K) being zero, is a guide to the eye only. 



Thus for large size system, the scaling law holds only within a very small parameter region. Nevertheless, our finding 
reveals the fact that there should exist a universal scaling law for such kind of phase transition. 

In the interaction-free case, from Section [3] we notice that the "PT symmetry phase hardly changes as the system 
size N increasing. In contrast, for medium interaction U, we investigate the phase diagram for finite size system 
by numerical simulation. The phase boundary is determined from the reality of the eigenvalues obtained by exact 
diagonalization. In Fig. [3] the phase diagrams are plotted as y c versus U for finite size chains. It shows that the 
on-site interaction breaks the VT symmetry drastically. Interestingly, there exist several !F*7~-symmetric regions and 
the number of such regions increases as N increases. The phase diagram shows rich structure in cases of medium U. 

4.4. Phase transition induced by bound-pair state 

In the following, we analytically investigate the boundaries of the !P7"~-symmetric phase of H. In strong on-site 
interaction case, there exists bound-pair band induced by the interaction U \ 45, 46[ 47]. In the presence of y, the VT 



symmetry is fragile. Here we focus on !PT-symmetric breaking phase transition caused by the bound-pair state. We 
analyse the phase diagram by introducing an effective Hamiltonian H e g, which describes the bound states of system 



H in strong on-site interaction U region. Based on the perturbation methods B47II . the effective Hamiltonian H e fi reads 
#eff = — (b]b M + h.c.) + I U + — I b]bi - — (b\bx + b' l N b N ) + 2iy (b\b } - b ] N b N ) , (54) 

where b^ = a^ 2 / V2 (b, - a 2 / V2) is the bound pair creation (annihilation) operator on site i. 

Similarly as in Section|4] the Hamiltonian H e s in single-particle invariant subspace can be solved by Bethe ansatz 
method. We are interested in the bound-pair band, which has the spectrum E (K) « [/+ (47 2 / U) (1 + cos K). Here, 
quasimomenta K of the bound-pair state satisfies the equation 

g(K) = 2(J 2 IUfsm(NK)(\ + cos K) + y 2 sin [(N - l)K] = 0. (55) 

As pointed above, together with the condition dg (K) /dK = 0, one can obtain the equation 

Z(K) = N(cosK + I) sin K -sin (NK) {cos (NK) + cos [(N -I) K]} = 0, (56) 
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Figure 5: The analytical boundaries i57\ , J5 8 1 . < I59 1 (red lines) and the phase diagram in large U region obtained by exact diagonalization for 
system H of size N = 10, 20, 40. The analytical boundaries are plotted in the region of U = 1 to 10, where / is set as the unit. It is noticed that the 
analytical boundaries fit the exact phase diagram well at the large U region. 



which solutions K c e (0, n) determine the boundaries of quantum phase for the effective Hamiltonian H s g. In other 
words, function Z (K) with N c zeros indicates [N c /2] + 1 unbroken !PT-symmetric regions of both H e g and H, where 
LA/ e /2J denotes the integer part of N c /2. To demonstrate this point, function Z (K) is plotted for different N in Fig. [4] 
One can see that, the number of solutions N c increases as N increases, which corresponds to the increasing unbroken 
VT -symmetric regions. To be specific, for N = 10, 20, and 40, it shows that N c = 1, 3, and 5, respectively. This 
indicates there are 1, 2, and 3 unbroken regions, which accords to the phase diagram in Fig. [5] Quantitatively, one can 
obtain the solutions of K c and y c from ( f55l l for these three cases numerically, which are listed in the following, 



K c /n = 
lny c /ln(J 2 /U) = 

K c /n = 
\ny c l\n(J 2 IU) = 

K c /n = 
lny c /\n(J 2 /U) = 



0.917579, 
1.515811, 

0.956567, 
2.124320, 

0.977429, 
2.746955, 



(for N = 10), 

0.938854, 
1.510840, 

0.971547, 
2.305626, 



0.914870, 
1.395938, 

0.955142, 
2.041551, 



(for AT = 20), 

0.942098, 0.933982, 
1.633456, 1.611533. 



(57) 
(58) 

(for N = 40). (59) 



According to the above analysis, equations (157115811591 represent the phase boundaries, which are plotted as red lines 
in Fig. [5]as comparison. It shows that the boundaries obtained from equations d57l l58l l59l are in agreement with the 
exact phase diagram well in large U regime. 

For large Af, the number of unbroken regions can be estimated as the integer around y/2N /n, the analytical results 
for different N of Fig. |4]is listed in ( f60b as comparison, it is noticed that the analysis accords with the plots in Fig. [4] 



N 10 
V2A77T 1.42 
[ V2A77T] 1 



20 
2.01 
2 



40 
2.85 
3 



200 
6.37 
6 



1000 
14.24 
14 



5000 
31.83 
32 



(60) 



flfN, correspondingly, the ceiling 



The ceiling (highest) phase boundary is determined by the root of d56l l near Kq 
phase boundary y c is near yo * (J 2 /U) y/2/ (N - 1). 

Now we turn to analyse the floor (lowest) phase boundary for H e ff , which corresponds to the solutions K c being in 
the region [(Af - l)n/N, n]. Setting K - n - {\ - 5) n/N with 6 <k 1 for large N, equation (TS6b can be reduced to 

(61) 



#{1 - cos [(1 - 6) n/N]} sin [(1 -6) n/N] - sin (nd) {cos (nS) - cos [nd + (1 - 6) n/N]} = 0. 
Moreover, by applying the Taylor expansion, it becomes 

2(N- l)6 2 + 36- 1 =0. 

Solving this equation, one can obtain 6 C = ( V8A^ + 1 - 3)/ (4Af - 4), and then the exceptional point y c as 

J 2 




Sc{s c -iy 

1)<5 C + 1 



(62) 



(63) 



We plot the exact and analytical approximation results for the floor phase boundary of H e ff, the exact result for the 
original Hamiltonian H in Fig. [6] as a comparison. The plots fit well, especially at large U case. Remarkably, it also 
gives a good approximation for those of medium U. 




Figure 6: Plots of the floor phase boundary of the phase diagram for N = 10,20,40. Black solid lines are the exact results obtained from 
diagonalization of the Bose-Hubbard Hamiltonian H. Red squares represent the exact result obtained from the effective Hamiltonian // c ff ■ Green 
circles are the analytical approximation ( 1631 . It is noticed that H e g and the approximation effectively describe the floor phase boundary of H. 



It is observed from Fig. |6]that // e tf gives quite good description of the phases of H for system with U > 1. On the 
other hand, the critical energy E c can be approximately expressed as 

E C «U+ T (\-6 c f (64) 

UN 2 



for large N. From the expression (163) , (1641 of y c , E c , we have 

J 2 E c -U \ 2 2Uy c 

l -w—)*^ (65) 

which is a universal scaling law for such a phase transition in large N, U case. Numerical simulations are performed 
to investigate the scaling behavior. We compute the quantities y c and E c for finite systems H e g , which are plotted in 
Fig. [7] as comparison with the analytical results (1631) . (l64l and ( |651 ). It shows for large N, they are in agreement with 
each other. 

In the limit N goes to infinity, we notice that y c = 1 when U = 0, means that in the region y < 1, system H 
is in the unbroken !P7~-symmetric phase. In non-zero interaction ( U + 0) case, from the above discussions we have 
the floor phase boundary y c « fin (UN) and ceiling phase boundary yo « ( J 2 /U) V2/(W - 1) for large N. As 




Figure 7: Plots of I63K <64t . {65) and the corresponding numerical simulation obtained by exact diagonalization of H e g . J = 1 is the unit. In (a), 
(b), the blue squares, circles and crosses indicate the numerical results for the cases of U = 1, 10 and 100, respectively. In (c), the red squares, 
green circles and blue crosses stand for U = 1, 10 and 100. The black lines are the plots of the corresponding analytical expressions j63h j64h 
)651 . It shows that they are in agreement in large N region. 
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goes to infinity, it is noticed that lim^M (y c ) — and lim^oo (yo) = for non-zero finite U. In other words, there 
always exist bound-pair bound states with conjugate complex energies, which result in the 'PT-symmetry breaking 
for non-zero y. 

System features experience huge changes as system parameter approaching the VT -symmetric phase transition 
point (exceptional point). Distinguished from phase transition in Hermitian system, the non-analytic properties in 
VT -symmetric phase transition are caused by the Hamiltonian becoming a Jordan block operator at the exceptional 
point. Correspondingly, dynamical behavior near the exceptional point experiences dramatic changes, e.g. the power 
oscillation amplitude becomes significantly large and the corresponding oscillation frequency becomes small. These 
features could be useful for weak signal detection, or for amplifier design. Phase diagram and the scaling behavior 
show us rich information where phase transition happens, this maybe helpful for the design and application of quantum 
devices. 



5. Conclusion and discussion 

In this paper, we have studied the scaling behavior and phase diagram of a non-Hermitian PT -symmetric Bose- 
Hubbard model. For interaction-free case, the metric operator is constructed, which is employed to investigate the 
particle number operator and the corresponding Hermitian counterpart. The derived properties of the metric operator, 
similarity matrix and equivalent Hamiltonian reflect the fact that they have a common feature: all the matrix elements 
change dramatically with diverging derivatives near the exceptional point. For nonzero U case, it is found that even 
small on-site interaction can break the VT symmetry drastically. It has been demonstrated that the scaling behavior 
can be established for the exceptional point in both small and large U limit. Based on numerical approach, we find 
that the phase diagram shows rich structure for medium U and analyse the phase transition boundary. Finally, it 
is worthwhile to point out that the phase transition discussed differs from the original quantum phase transition in 
a Hermitian system fl49ll . The former aims at certain eigenstates of a non-Hermitian Hamiltonian, while the later 
concerns only the ground state of a Hermitian one. However, our finding reveals that both of them exhibit the scaling 
behavior, which may be due to they both relate to the spontaneous symmetry breaking. 

Finally, we would like to discuss the relevance of the present model to a real physical system. The Hermitian 
Bose-Hubbard model is the simplest model capturing the main physics of not only cold atoms in optical lattice also 
photons in nonlinear waveg uide 00 SB. The effective non-Hermitian Bose Hubbard Hamiltonian //loss, is 
introduced when the closed system couples to the continuum [5(i 51, 52, 53]. A experimental realization of such 
an open system could be achieved by tunneling escape of atoms from a magneto-optical trap [41] or using lossy 
cavities 15411. The Hamiltonian //loss has been investigated by solving the master equation 14 lll55[ 15611 or Schrodinger 
equation 137 , 38]. Although the present model Hamiltonian H (Q} contains an extra gain term iy, it has connection to 



a Hamiltonian //loss by applying a constant energy shift: 



H 



Loss 



H 



N 



(66) 



The dynamics of //loss can be obtained by solving the following master equation 



N-l 

= -i [H ,p] - y^(a!fl/p + pa\cii - 2a;pa\)- 2y{a ] N a N p + pa ] N a N - 2a N pa} N ), 

1=2 



(67) 



where p is the density matrix of //loss, Ho = -J J] 1 ^ 1 (ajai+i+H.c.) + (U/2) aV'aJ. 

Based on the result of this paper, some important features of //loss can be observed. Firstly, although Hamiltonian 
//loss is not invariant under PT transformation, the eigenstates of Z/l OS s are st iU fT symmetric within the unbroken 
region of H. Secondly, two Hamiltonians H and Z/l OS s share the same dynamics except the extra decaying factor, i.e., 
t//(t) — > il/(t)e~ 7 '. From this perspective, the phase boundary as well as the scaling law presented in this paper, can 
be observed from the dynamics of //loss- As a future work, it is interesting to compare the results obtained by two 
methods. Actually, it has been explored for a two-site example |55, 56ll . 
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Appendix A. Operators hi and h kl 

This appendix provides the examples to demonstrate that operators hi = crfai and h k , — a[a.k, are not observables. 
We consider the single-particle case as an illustrative example, the metric operator can be expressed as 



»/ = Z «*(**)* WW. 

and the operators are hi = \l) (l\ and h kl = \k/) (ki\. Then we have 

{qhi - h\ j]) hh = (ji | {rjhi - h\rf) \j 2 ) 



(A.1) 



(A.2) 



or more explicitly 



(j]hi-h\ri) hn 



-iz k g' k (g j k 2 )\ <ji = i±& 

h=l = h 

or h 



(A.3) 



0, 



We note that ^(guYgf do not vanish in general case, which can be seen in the following illustrative example. We 
consider the case with N = 2, from the operator 772 listed in the Table [T] we notice that 



mn\ - Km = 



m n 2 - n 2 T] 2 



which means mh[T] 2 ' ^ "1 ■ Accordingly, it leads to 



i 
yy 2 -y 2 \ i 
i 



i 



*0, 



(A.4) 



T]{lhi)r] 1 = IrjhiT] 1 + Ihirji] 1 = (/«;)', 



(A.5) 



which shows that the position operator lh\ is not an observable. This accords with the conclusion of Ref. 
Similarly, the operator ?\, for N - 2 system in single-particle case, has the form 



i^r—f- 
1 



-./ 



in coordinate space. Straightforward algebra shows 

mhk x - «], m 

mn kl - hi m 



■^J 2 - y 2 - iy ^ 
■^J 2 —y 2 + iy -J 

J yjj 2 -y 2 + iy 

^J 2 - y 2 - iy J 



(A.6) 



-i 

y[j^f-\ i 
i 



*0, 



(A.7) 



yp^\ -i 

which means mh^rf^ + ft k . Then we conclude that operator is not an observable. 
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